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1  Introduction 


We  report  here  on  the  work  performed  during  the  three  years  (April  1st,  2012  -  March  31st, 
2015)  of  the  contract  N00014-12-1-0383  on  Multiscale  materials  science:  a  mathematical  ap¬ 
proach  to  the  role  of  defects  and  uncertainty. 

The  presence  of  numerous  length-scales  in  material  science  problems  represents  a  daunt¬ 
ing  clialleiige  for  numerical  simulation.  Quantifying  the  effects  of  defects,  and  more  generally 
any  uncertainty  arising  from  data,  discretization,  and  the  mechanical  model  for  an  associated 
numerical  method  has  become  an  increasingly  important  aspect  of  multiscale  analysis.  Uncer¬ 
tainty  in  materials  science  problems  assumes  various  forms,  and  includes  defects  in  crystals, 
impurities  or  heterogeneities  in  continuous  media,  . . .  Our  aim  is  to  develop  new  mathemati¬ 
cal  and  numerical  tools,  including  probabilistic  approaches,  to  address  the  current  challenging 
problems  of  interest  in  materials  science. 

Ill  this  report,  we  first  focus  on  diiveloppiiig  affordable  mmierical  methods  in  the  context 
of  stochastic  homogenization. 

Many  partial  diffiuential  eciiiatioiis  of  materials  seii'iiee  indeed  involve  highly  oscillatory 
coefficients  and  small  kuigth-scales.  Homogeiiizatioii  thixiry  is  concerned  with  the  derivation 
of  averaged  equations  from  the  original  oscillatory  equations,  and  their  treatment  by  adequate 
numerical  approaches.  Stationary  ergodic  random  problems  (and  the  associated  stochastic 
homogenization  theory)  are  one  instance  for  modelling  uncertainty  in  continuous  media.  The 
theoretical  aspects  of  these  problems  are  now  well-understood,  at  least  for  a  large  variety  of 
situations.  On  the  other  hand,  the  numerical  aspects  have  received  less  attention  from  the 
mathematics  community.  Standard  methods  available  in  the  literature  often  lead  to  very,  and 
soiiK'tiiiios  prohibitively,  costly  computations,  while  affordable  approaches  do  not  rest  upon 
a  sound  mathematical  theoretical  setting. 

Our  aim  is  to  compute  the  effective  or  houiogenizc'd  coefficic'iits  that  represent  the  be¬ 
havior  of  the  material  at  a  macroscopic  scale.  In  a  first  part  (see  Section  2),  we  give  a  very 
brief  overview  of  classical  results  of  stochastic  homogenization  theory.  Then  we  explain  why 
such  results  lead  to  practical  difficulties.  In  particular,  though  the  effective  properties  of 
the  material  considered  here  are  deterministic,  their  approximations  require  solving  partial 
differential  eriuatioiis  with  random  coefficients  posed  on  large  domains.  This  induces  very 
expensive  computations. 

In  Section  3,  we  consider  one  possible  technique  to  decrease  the  cost  of  such  computations, 
namely  the  control  variate  technique.  This  variance  reduction  technique,  which  is  based 
on  using  a  surrogate  model,  is  classical,  and  has  already  been  successfully  applied  in  other 
scientific  fields.  We  show  here  that  it  can  also  be  used  in  the  general  framework  of  stochastic 
homogenization.  In  our  case,  the  surrogate  model  that  we  use  is  inspired  by  a  defect-type 
theory,  where  a  perfect  periodic  material  is  perturbed  by  rare  defects.  This  model  has  been 
introduced  in  [9]  in  the  context  of  weakly  random  models.  Here,  we  address  the  fully  random 
case,  and  show  that  the  perturbative  approach  proposed  in  [9,  11]  can  be  turned  into  an 
efficient  control  variable.  When  the  perturbation  is  not  small,  the  perturbative  approach, 
although  not  accurate,  is  good  enough  to  reduce  the  variance  of  a  complete  computation.  For 
a  large  class  of  initial  highly  oscillatory  coefficients,  the  resulting  approach  yields  a  significantly 
moix'  precise  estimation  of  the  properties  of  the  effective  material,  at  equal  computational  cost, 
than  Monte  Carlo  approaches  or  other  variance  reduction  approaches  that  we  have  previously 
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adapted  to  the  homogenization  setting. 

In  Section  4,  we  address  questions  related  to  non-periodic  modelling  of  multiseale  ma¬ 
terials.  Random  modelling  is  indeed  a  standard  modelling  response  to  situations  when  the 
idealized,  say  periodic,  modelling  is  inappropriate.  It  is  however  commonly  admitted,  and 
observed  in  practice,  that  random  modelling  leads  to  possibly  prohibitively  computation¬ 
ally  expensive  problems.  Then^  is  a  definite  theoretical  and  practical  interest  in  generalizing 
modelling  of  perfect  materials  in  directions  different  from  the  random  paradigm. 

Within  the  current  contract,  we  have  investigated  this  idea,  which  has  not  been  explored 
in  the  previous  ONR  grant.  The  general  approach,  that  we  describe  in  Section  4,  consists  in 
approximating  at  the  fine  scale  the  solution  to  an  elliptic  eciuatioii  with  oscillatory  coefficient 
when  this  coefficient  consists  of  a  “nice”  function  (say  periodic)  which  is  perturbed  by  a  local 
defect.  To  date,  the  questions  that  have  been  investigated  are  essentially  of  a  theoretical 
nature.  They  yet  give  rise  to  an  associated  numerical  endeavour,  that  is  now  on  the  horizon. 

In  Section  5,  we  also  investigate  a  research  direction  that  has  not  been  explored  in  the 
previous  ONR  grant.  All  models  that  involve  a  random  parameter  require  (at  least  a  partial, 
and  most  often  a  complete)  knowledge  of  the  distribution  of  this  random  parameter.  Now, 
access  to  this  distribution  is  practically  difficult.  One  is  therefore  bound  to  assume  a  given 
form  (Gaussian,  . . . )  for  the  distribution  and  proceed  with  the  computation.  A  question  of 
major  practical  interest  is  to  a  posteriori  prove,  or  disprove  the  validity  of  this  assumption. 
Put  differently,  tests  of  hypotlieses  in  the  c  ontext  of  engineering  problems  is  an  important 
issue.  A  preliminary  step,  before  trying  to  identify  the  distribution  of  the  random  parameters, 
is  to  assume  a  specific  form  of  this  distribution,  dc'peiidiiig  on  a  few  (piantities  (e.g.  assume  a 
Gaussian  distribution  with  an  unknown  variance)  and  identify  these  quantities.  The  parameter 
identification  problem  that  we  consider  in  Section  5  exactly  fits  in  this  preliminary  step. 

The  works  described  below  have  been  performed  (jointly  or  separately)  by  Claude  Le  Bris 
(PI),  ErARrie  Legoll  (Co-PI)  and  William  Miiivielle  (Ph.D.  student). 


2  Basics  of  stochastic  homogenization 

[A  detailed  presentation  can  be  read  in  [8,  14]-] 

For  the  sake  of  completeness,  and  to  fix  notations,  we  give  here  a  very  brief  overview 
of  classical  results  of  stochastic  homogenization  theory.  The  reader  familiar  with  stochastic 
homogenization  can  proceed  directly  to  our  contributions,  detailed  in  Sections  3,  4  and  5. 

Stochastic  homogenization  is  best  understood  in  the  light  of  the  easiest  context  of  ho¬ 
mogenization:  periodic  homogenization.  This  is  the  reason  why  we  begin  by  laying  some 
groundwork  in  the  periodic  context,  before  turning  to  stochastic  homogenization  per  se  in 
Section  2.2.  We  refer  to,  e.g.,  the  monographs  [17,  19,  22]  for  more  details  on  homogenization 
theory,  and  to  the  review  article  [8],  where  we  addressed  some  computational  challenges  in 
numerical  stochastic  homogenization.  A  super  elementary  introduction  is  contained  in  [15]. 
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2.1  Periodic  homogenization 

For  consistency,  we  recall  here  some  basic  ingredients  of  periodic  homogenization  theory.  We 
consider,  in  a  regular  bounded  domain  V  of  the  problem 


^per 


=  /  in  V, 


(1) 


where  the  matrix  Aper  is  symmetric  definite  positive  and  Z'^-periodic.  We  manipulate  for 
simplicity  symmetric  matrices,  but  the  discussion  carries  over  to  non  symmetric  matrices  up 
to  slight  modifications. 

The  microscopic  problem  associated  to  (1),  called  the  corrector  problem  in  the  terminology 
of  homogenization  theory,  reads,  for  p  fixed  in 


-div  [Aper(y)  {p  +  Vwp{y))]  =0  in 
Wp  is  Z'^-periodic. 


(2) 


It  has  a  unique  solution  up  to  the  addition  of  a  constant.  Then,  the  homogenized  matrix  A* 
is  such  that,  for  any  p  € 


A*p=  /  Aper(y)  (p  + Vu;p(y))(iy, 

■Jq 

where  Q  =  (0,  is  the  unit  cube.  The  main  result  of  periodic  homogenization  theory  is 
that,  as  e  goes  to  zero,  the  solution  to  (1)  converges  to  u*  solution  to 

(  —div  =  /  in  V, 

\  (3) 

n*  =  0  on  dT>. 

Several  other  convergences  on  various  products  involving  Aper  and  also  hold.  All  this 
is  well  documented. 

The  practical  interest  of  the  approach  is  evident.  No  small  scale  e  is  present  in  the 
homogenized  problem  (3).  At  the  price  of  only  computing  d  periodic  problems  (2)  (as  many 
problems  as  dimensions  in  the  ambient  space),  the  solution  to  problem  (1)  can  be  efficiently 
approached  for  e  small.  A  direct  attack  of  problem  (1)  would  require  taking  a  meshsize 
smaller  than  e.  The  difficulty  has  been  circumvented.  Of  course,  many  improvemeiits  and 
alternatives  exist  in  the  literature. 

2.2  Stochastic  homogenization 

Stochastic  homogenization  has  a  mathematical  setting  that  is  more  involved  than  that  of  the 
periodic  case. 

We  consider  the  theoretical  setting  of  stationary  ergodic  homogenization,  with  a  discrete 
shift  operator,  which  intuitively  means  the  following.  Pick  two  points  x  and  y  ^  x  at  the 
microscale  in  the  material  and  assume  y  =  x-t-k  with  k  G  Z'^.  The  particular  local  environment 
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seen  from  x  (that  is,  the  microstructure  present  at  x)  is  g(‘iieric-ally  different  from  what  is  seen 
from  y  (that  is,  the  microstructure  present  at  y).  However,  the  average  local  environment 
in  X  is  assumed  to  be  identical  to  that  in  y  (considering  the  various  realizations  of  the  random 
material).  In  mathematical  terms,  the  probability  law  of  the  microstructures  is  the  same. 
This  is  stationarity.  On  the  other  hand,  ergodicity  means  that  considering  all  the  points  in 
the  material  is  eciuivaleiit  to  fixing  a  point  x  in  this  material  and  considering  all  the  possible 
microstructures  present  there. 

Consider  now  the  problem 


I  _div(.4(^  C.)V»')=/  ta  ®, 

u®  =  0  on  dV, 

where  the  matrix  encoding  the  properties  of  the  material,  is  assumed  to  be  stationary  (or, 
equivalently,  statistically  homogeneous).  The  solution  u®  to  (4)  converges,  when  £  — >■  0,  to 
the  solution  u*  to  (3)  where  the  homogenized  matrix  is  now  given,  for  any  p  G  by 


=  E  /  A{y,-){pAVwp{y,-))  dy 


Q 


and  the  corrector  problem  now  reads 

-diy[A{y,uj){p  +  VWp{y,u}))]=0  in 


Vwp  is  stationary,  E 


'Q 


Vwp{y,-)dy'^  =  0. 


(5) 


(6) 


A  striking  differenc-e  between  the  stochastic  setting  and  the  p(^riodic  setting  can  be  observed 
comparing  (2)  and  (6).  In  the  periodic  setting,  the  corrector  problem  is  posed  on  a  bounded 
domain,  namely  the  periodic  cell  Q.  In  sharp  contrast,  the  corrector  problem  (6)  of  the 
random  setting  is  posed  on  the  whole  space  and  cannot  be  reduced  to  a  problem  posed 
on  a  bounded  domain. 

The  fact  that  the  random  corrector  problem  is  posed  on  the  entire  space  has  far  reaching 
consequences  for  numerical  practice.  In  particular,  this  is  the  reason  why  standard  methods 
available  in  the  literature  often  lead  to  very  costly  computations,  a  fact  which  in  turn  motivates 
our  work. 


3  Variance  reduction  techniques  for  stochastic  homog¬ 
enization 

[Work  expanded  in  [6].] 

3.1  The  direct  numerical  approach 

Practical  approximations  of  the  homogenized  matrix  in  random  homogenization  are  not  easily 
obtained,  owing  to  the  fact  that  the  corrector  problem  (6)  is  set  on  the  entire  space.  In 
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practice,  truncations  have  to  be  considered,  and  the  exact  homogenized  coefficients  are  only 
obtained  in  the  asymptotic  regime. 

In  practice,  the  matrix  A*  is  indeed  approximated  by  deffined  by 

A%{u)p=-^f  A{y,u)  {p+VWp{y,uj))  dy,  (7) 

IVVI  ./Qjv 

which  is  obtained  by  solving  the  corrector  problem  on  a  truncated  domain,  say  the  cube 
=  (-iV,  nY  C 


I  -div  (p+ Vn;^(-,w))]  =  0  on  E^, 

wY(-,oj)  is  QAT-periodic. 

Although  A*  itself  is  a  deterministic  object,  its  practical  approximation  is  random,  for 
all  N  finite.  It  is  only  in  the  limit  of  infinitely  large  domains  Qn  that  the  deterministic  value 
is  attained.  It  has  indeed  been  shown  in  [18,  Theorem  1]  that  limAr_,.oo  =  A*. 

We  now  remark  that  the  error  can  be  decomposed  as 

-  /I*  =  -  E  (>!*„))  +  (E  (>l*„)  -  A*) , 

where  the  first  term  is  a  statistical  error  and  the  sec-ond  term  is  a  bias  (systematic  error).  We 
focus  here  on  the  statistical  error,  and  propose  approaches  yielding  better  approximations  of 
E[A)^],  for  a  given  truncated  domain  Qn-  Optimal  estimates  on  the  variance  of  have 
been  established  in  [23,  Theorem  1.3  and  Proposition  1.4].  See  also  [20,  Theorem  1].  In  these 
works,  it  has  been  noted  that  the  systematic  error  is  much  smaller  than  the  statistical  error,  in 
the  sense  that  the  latter  decays  with  a  slower  rate  with  respect  to  N  than  the  former.  These 
results  are  consistent  with  numerical  observations.  For  large  values  of  N,  the  statistical  error 
(that  we  address  here)  is  therefore  dominating  over  the  systematic  error.  There  is  thus  a 
dehnite  practical  interest  to  focus  on  that  error  and  design  approaches  to  reduce  it,  as  we  do 
here. 

To  approximate  E[A^],  a  standard  way  is  the  Monte  Carlo  method.  We  give  ourselves 
a  set  of  M  independent  copies  (or  realizations)  (A"*)j<^<^^  of  the  random  coc-fficient  A. 
The  corresponding  truncated  problems  (8)  are  solved,  which  provides  us  with  a  sequence 
of  independent  and  identically  distributed  homogenized  matrices  A*Y^{ui),  defined,  for  any 
p  G  E'^,  by 

^TMp  =  A-t  f  (p+ i^)). 

IVAfI  Jqn 

where  w^’’^  is  the  solution  to  the  corrector  problem  (8)  associated  to  A"*.  Tluui  we  dehne 
the  empirical  mean 

1 

d-M  (An)  =  a*/" .  (9) 

m=l 

Since  the  matrices  A*^  are  i.i.d.,  the  strong  law  of  large  numbers  applies: 

/^A/  {A%)  (^)  — ^  E  (A^)  almost  surely. 

M— >-+oo 
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The  central  limit  theorem  then  yields,  for  any  component  1  <  i,  j  <  d, 

'/M  (pm  -  E  ((/i;,|«))  j  Ar(0, 1),  (10) 

where  the  convergence  holds  in  law,  and  jV(0, 1)  denotes  the  standard  Gaussian  law.  From  (10) 
it  is  commonly  admitted  that  the  exact  mean  E  ( j  lies  in  the  coiifidcuice  interval 


Jvav 

l‘M  ([.4;,|y)  -  l,96.i - ^ -  ,  PM  ([/lj,]y)  +  1.96 

The  value  Hm  thus,  for  both  M  and  N  sufficicnitly  large,  oft('n  adopted  as  an 

approximation  of  the  exact  value  [^*]y-  The  overall  computation  described  above  is  thus  very 
expensive,  because  each  realization  requires  a  new  solution  to  the  problem  (8)  of  presumably 
large  a  size  since  N  is  taken  large. 

In  the  sequel,  we  show  that,  using  a  control  variate  approach,  we  can  design  a  practical 
approach  that,  for  any  finite  N,  allows  to  compute  a  better  approximation  of  E  than  (9). 
Otherwise  stated,  for  an  equal  computational  cost,  we  obtain  a  more  accurate  (i.e.  with  a 
smaller  eonhdc'iice  interval)  approximation. 

We  detail  in  Section  3.2  the  general  principle  of  the  approach,  and  show  how  to  apply  it 
to  the  homogenization  setting  in  Sections  3.3  and  3.4.  We  next  collect  our  numerical  results 
in  Section  3.5  and  our  theoretical  results  in  Section  3.6. 


3.2  The  control  variate  technique:  general  principle 


Befor(5  applying  the  approach  to  our  specific  setting,  we  bri('fly  describe  here  the  control  variate 
approach  in  a  general  context.  Considering  a  random  variable  X,  our  aim  is  to  compute  its 
expectation  E(X).  In  the  sequel,  we  will  use  that  approach  for  the  random  variable 
for  any  entry  1  <  i,  j  <  d. 

As  point('d  out  above,  a  first  possibility  is  to  resort  to  M  i.i.d.  realizations  of  X,  denoted 
X’^(uj)  for  1  <  m  <  M.  The  expectation  is  then  approximated  by  the  Monte  Carlo  empirical 
mean 

m=l 

and  we  know  that,  with  a  probability  equal  to  95  %,  E  [X]  lies  in  the  conhdence  interval 


rMC 

-'a/ 


1.96 


\fM 


rMC 

-'a/ 


+  1.96 


Vm 


(11) 


To  reduce  the  variance  of  the  estimation,  consider  now  a  random  variable  Y,  the  expectation 
of  which  is  analytically  known.  Then,  for  any  deterministic  parameter  p  to  be  hx('d  later,  wc^ 
consider  the  controlled  variable 

D„(ui)  =  x(uj)  -  p(r(w)  -  E[y|) .  (12) 
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Since  E[F]  is  known  exactly,  sampling  realizations  of  Dp  amounts  to  sampling  realizations 
of  X  and  Y .  We  obviously  have  E[Dp]  =  E[X].  To  approximate  E[X],  the  control  variate 
approach  consists  in  performing  a  standard  Monte  Carlo  approximation  on  Dp.  We  hence 
consider  M  i.i.d.  realizations  of  Dp,  denoted  D'^{u)),  introduce  the  empirical  mean 

M  ..  M 

w = 17  e  ’ 

m=l  m=l 

and  write  that,  with  a  probability  equal  to  95  %,  E[Dp]  =  E  [X]  lies  in  the  confidence  interval 


/OV  _ 

y/M 


I?7  +  1.96 


■5 


\/M 


(13) 


If  p  and  Y  are  such  that  Var  [Dp]  <  Var  [X],  then  the'  width  of  the  above'  c'onfideiice  interval 
is  smaller  than  that  of  (11),  and  hence  we  have  built  a  more  accurate  approximation  of  E  [X]. 

The  choice  of  Y  is  problem-dep('iid('iit  and  is  discussed  Ix'low  in  our  speeific  context  (see 
Section  3.4).  Assuming  for  now  that  Y  is  given,  we  detail  here  the  classical  rationale  and 
technique  employed  to  choose  p  in  (12).  We  wish  to  pick  p  such  that  the  variance  of  Dp  is 
minimal.  Writing  that 


Var[Dp]  =  Var[X]  -  2pCov[X,  T]  +  p^VarfT], 
we  see  that  the  optimal  value  of  p  reads 

,  ,,  ,  ,  Cov[x,y] 

p  =  argmm  Var[Dp]  =  . 

For  this  choice,  we  have,  using  the  Cauchy-Schwarz  inequality, 

(Cov[x,y])2 


(14) 


Var  [Dp*]  =  Var[X]  1 


Var  [X]  Var  [T] 


<  Var[X]. 


We  thus  observe  that,  for  any  choice  of  Y ,  we  can  choose  p  such  that  the  variance  of  Dp 


is  indeed  smaller  than  that  of  X.  Of  course,  the  ratio  of  variances 


Var  [Dp.. 


-,  which  is  di- 


Var[X] 

rectly  related  to  the  gain  in  accuracy,  depends  on  Y ,  and  more  precisely  on  the  value  of 

— —  ’ — r^.  The  larger  the  correlation  between  X  and  Y ,  the  better.  In  particular,  we 

Var[X]Var[y  ] 

see  that  the  control  variable  Y  needs  to  be  random. 

In  practice,  we  do  not  have  access  to  the  optimal  value  (14),  which  involves  exact  expec¬ 
tations.  One  possibility  (which  is  the  one  we  adopt  here)  is  to  replace  (14)  by  the  empirical 
estimator 
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1  ^ 

where  IjLxi{X)  =  —  X"^{u)).  This  choice  corresponds  to  minimizing  with  respect  to  p  the 


m=l 


1  ^  2 

empirical  variance  of  defined  as  —  ^  -  Pm{X))\  where  D^{uj)  =  X-(a;)  - 


m=l 


p(y'^{u))  —  E[T]j.  The  expectation  E(JT)  is  then  approximated  by 


1 

=  mJ2  ■ 

m=l 

3.3  A  weakly  random  setting:  rare  defects  in  a  periodic  structure 

As  pointed  out  in  the  introduction,  the  surrogate  model  that  we  use  to  build  our  controlled 
variable  is  inspired  by  a  defect- type  model,  introduced  in  [9,  10,  11]  in  the  context  of  weakly 
random  models,  and  which  we  describe  now. 

3.3.1  Presentation  of  the  model 

Assume  that,  in  (4),  the  random  matrix  A  is  of  the  form 

A(x,  Cu]  Aj^{x^ui^  Ap0i' 4“  67^  Cu)  ^C^per  (^)  Apei'(x)^  (fb) 

where  Aper  and  Cper  are  ZAperiodic  matrices,  and 

^(a;,a;)  =  ^  lQ+k{x)Bl{uj),  (17) 

fcez*^ 

where  are  i.i.d.  scalar  random  variables.  We  furthermore  assume  that  B^  follows  a 

Bernoulli  law  of  parameter  ?]  €  (0, 1): 

V(Bl  =  1)  =  r,,  V(Bl  =  0)  =  1  -  r,.  (18) 

In  each  cell  Q-\-k,  the  field  A  is  equal  to  Aper  with  the  probability  1  —  ?],  and  equal  to  Cper  with 
the  probability  p.  When  p  is  small,  then  (16)-(17)-(18)  models  a  periodic  material  (described 
by  Aper)  that  is  randomly  perturbed  (and  then  described  by  Cper).  The  perturbation  is  rare 
when  p  is  small  (therefore  the  material  is  described  by  Aper  “most  of  the  time”),  and  thus  it 
can  be  considered  as  a  defect.  However,  the  perturbation  Cper  —  Aper  is  not  small.  We  refer 
to  [11]  for  practical  examples  motivating  this  framework. 

On  Fig.  1,  we  show  two  realizations  of  the  field  Arf{x,oj)  (on  the  domain  for  N  =  20) 
for  soiiK^  specific  choices  of  Aper  and  Cper  (see  [11,  Fig.  4.2]  for  more  details).  On  the  right 
part  of  that  hgure,  we  s<!t  p  =  0.4,  which  is  close  to  the  value  p  =  1/2,  when  defects  are  as 
frequent  as  non-defects. 

Note  that  specifying  Ajj{x,  co)  on  Qn  simply  amounts  to  specifying  the  values  of  B^iuj)  for 
all  k  such  that  k  +  Q  C  Qn- 
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Figure  1;  Two  instances  of  material  (16).  Left  (ry  =  0):  perfect  material  with  circular 
inclusions  located  on  a  periodic  network.  Right  (ry  =  0.4):  perturbed  material  (each  inclusion 
is  deleted  with  a  probability  equal  to  0.4). 

3.3.2  Weakly-random  homogenization  result 

Consider  the  model  (16)-  (17)-  (18).  The  random  variable  can  take  only  two  values,  0 

or  1.  Therefore,  on  the  domain  Qat,  tlic're  are  only  a  finite  number  of  rcializations  of  Aj,{x,u). 
The  realizations  with  the  highest  probabilities  are  as  follows. 

With  probability  (1  —  there  are  no  defects  in  Qv,  and  the  realization  actually 

corresponds  to  the  perfect  periodic  situation.  We  introduce  the  periodic  corrector  Wp,  solution 
to 

-div  (^per  {p  +  Vu’p))  =  0,  Wp  is  Q-periodic,  (19) 

and  the  associated  matrix  ^pgj.,  obtained  by  periodic  homogenization: 

Vp  G  M",  a;,,  p  =  [  Tlper  (p  +  VrcJ)  .  (20) 

Jq 

With  probability  p(l  —  there  is  a  unique  defect  in  Qw,  located,  say,  in  the  cell  k  +  Q 

(sc-e  F'ig.  2).  Lc't  us  defiiu' 

=  ^per  +  Ifc+Q  ^Cper  —  ^per^  ,  (21) 

the  associated  corrector  solution  to 

—div  (p  +  Vwp’^’'^))  =  0,  Wp'^'^  is  Q^-periodic,  (22) 

and  the  homogenized  matrix  Ah^j^,  given  by 

VpeR",  4(p  +  V™p-«),  (23) 

With  probability  p^(l  —  there  are  two  defects  in  Qn,  located,  say,  in  the  cells  k  +  Q 

and  I  +  Q  (see  Fig.  2).  Let  us  define 

-d.2  .Tper  +  ^lfc-)-Q  +  1;-|_Q^  ^Cper  —  j4per^,  (24) 
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the  associated  corrector  solution  to 


-div  +  =0, 

and  the  homogenized  matrix  given  by 

VpGM'^,  Alk,i,NP^ 


w. 


2,k,l,N 


is  Q^vP^riodic, 


A’^'^  {p  +  . 


iQivI  Jqn 

All  I.Ik^  other  eonliguratioiis  (with  three  defects  or  more)  have  a  smaller  probability. 


(25) 


(26) 


•  ••••  •••• 
•  •  ••••••• 


Figure  2:  Left;  material  modelled  by  A'^,  with  a  single  defect.  Right:  material  modelled  by 
A2’'',  with  two  defects. 


As  formally  shown  in  [11,  Section  3.2],  we  then  have  the  following  result,  for  any  given  N: 


E  [Al^]  =  A;^,  +  r/Af  +  t^^A;  +  0{v^), 


7-rN 


where  A*^^  is  given  by  (20)  and 


-  X!  ^ 


•k 

per 


A 


=  5  E 


k^l^Xj^jk^l 


-^2,k,l,N  ^l,k,N  +  ^per^  i 


where 

We  note  that 


Xjv  : —  |/>;  G  Q  k  CL  Qn  | . 


-^N  -^k,N  ,  -rN  1  ^  -Tk,l,N 

A=2^Adefi  and  A2=-2^  2^  Asdef, 

/cGXjv  k€XN^^XN,l^k 


(27) 


(28) 


where  A^f^f  (resp.  A^dL^)  is  the  “marginal  contribution”  to  the  homogenized  matrix  from  a 
conhguration  with  a  single  defect  in  A:  +  Q  (resp.  two  defects  in  A;  +  Q  and  I  +  Q): 


/A'^  —  A*  _  A* 

^Idef  ^per’ 


-^kJ,N 


a: 


2  def 


~  ■^2.k,l,N  ^l,k,N  ^ 


per* 


(29) 

(30) 
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We  refer  to  [11]  for  illustrative  numerical  results. 

Due  to  periodic  boundary  conditions  (22),  that  are  reminiscent  of  the  periodic  boundary 
conditions  in  (8),  we  have  that 

does  not  depend  on  k.  (31) 

Likewise,  i^  ,  v  depends  only  on  k  —  1.  Thus,  there  is  only  one  problem  (22)  to  be  solved 
(say  for  k  =  Q).  Likewise,  there  are  |I/v|  —  1  problems  (25)  to  be  solved  (say  for  k  =  0  and 
1^0),  and  not  \Xn\  {\Xn\  ~  !)•  Noticing  that  (25)  is  a  problem  parameterized  by  I,  the  authors 
of  [16]  have  shown  how  to  use  a  Reduced  Basis  approach  to  further  speed-up  the  computation 

— JV  .  .  — N 

of  A.2  .  In  practice,  one  can  still  obtain  a  good  approximation  of  A2  without  solving  all  the 
\Tn\  —  1  problems  (25). 


3.4  Application  to  stochastic  homogenization 

We  now  introduce,  for  the  model  (16)-(17)-(18),  a  control  variate  approach.  Our  aim  is 
now  to  address  the  regime  when  rj  is  not  close  to  0  or  1  (the  approximation  (27)  is  therefore 
not  accurate  enough).  Recall  also  that,  in  view  of  the  discussion  of  Section  3.2,  we  need  a 
random  surrogate  model  to  build  our  controlled  variable.  In  what  follows,  we  first  build  an 
api)roxiinate  model  based  on  configurations  with  a  single  defect  (see  Section  3.4.1),  and  next 
turn  to  building  a  better  approximate  model  that  also  uses  configurations  with  two  defects 
(see  Section  3.4.2).  As  will  be  seen  below,  this  second  approximate  model  not  only  depends 
on  the  quantity  of  defects,  but  also  on  their  geometry,  that  is  on  where  the  defects  are  located 
in  Qjv- 

Numerical  results  obtained  with  these  approaches  are  collected  in  Section  3.5.  On  the 
particular  test-case  at  hand,  we  will  observe  (see  Fig.  5  below)  that  the  approach  using  the 
single  defects  provides  a  variance  reduction  ratio  close  to  6,  while  the  approach  using  single 
defects  and  pairs  of  defects  provides  a  gain  roughly  6  times  larger  (here  of  the  order  of  40). 


3.4.1  A  first-order  model 


Recall  (see  (16)-(17))  that 

A(x, cc)  A-q{x^uj^  Apei.(37)  -{- 6,^ (x, cc)  ^(17p0[. (37)  Ap0j.(x)^ 

where  Aper  and  Cper  are  Z'^-periodic  matrices,  and 


6„(x,w)  =  lQ+k{x)Bl{oj), 
k&Z'i 

where  are  i.i.d.  scalar  random  variables.  Introduce 

E  wAUf.  (32) 


— k,N 

where  d('hned  by  (29),  is  the  “marginal  contribution”  to  the  homog('uiz('d  matrix  coming 

from  tlu'  couhguration  with  a  single  defect  located  in  k  A  Q.  In  view  of  (28),  we  notice  that 


E 


A? 


N 


=  E  E  K|  A:Z,  =  riJ2^< 


‘k,N 

def 


r/Af, 


12 


which  is  the  first  order  correction  in  the  expansion  (27).  When  r]  is  small,  the  expectation  of 
+  is  a  good  approximation  of  the  expectation  of  accurate  up  to  an  error 

of  the  order  of  r]^.  Thus  +  (co)  is  a  good  surrogate  model  for  A*  j^{u).  Following  (12), 
we  now  introduce  our  controlled  variable  as 


-  f  ('^i'''(‘^)  -  W)  ■ 


(33) 


In  view  of  (32),  (29)  and  (31),  we  recast  (33)  as 


D];^{u)  =  Al^{u)-p 


n\^N\ 


^Idef- 


(34) 


Remark  1  Note  that,  in  (34),  A*j,^{<jo)  and  S^(cn)  are  correlated.  Indeed,  in  practice, 

we  start  by  drawing  a  realization  of  the  random  variables  for  all  k  E  In .  This  deter¬ 
mines  first  Efceijv  second  the  field  A{x,u)  on  Qn,  from  which  we  compute  the 

associated  A*  ,,^{u)  following  (7)  (8).  The  fact  that  A*  and  correlated 

is  important  in  view  of  variance  reduction,  as  explained  in  Section  3.2. 


Computing  M  realizations  of  D];^[u))  therefore  amounts  to: 

•  (offline  stage)  determine  by  solving  the  problem  (19)  (20)  on  Q  and  solving  only 

once  the  problem  (22)-(23)  on  Qn  (say  for  k  =  0),  and 

•  (online  stage)  solve  M  corrector  problems  (7)  -(8)  on  Qn  (for  M  i.i.d.  realizations  of  A 
on  Qn),  and  evaluate  Dj;^{u)  according  to  (34). 

Let  Cn  denote  the  cost  to  solve  a  single  corrector  problem  on  Qn-  The  Monte  Carlo  empirical 
('Stimator  and  the  Control  Variate  empirical  (estimator,  defiiuxi  respectively  by 


rMC 

•'Af 


m=l  m=l 


therefore  share  the  same  cost  (Af  Cat  for  the  former,  (1  +  M)  Cpf  for  the  latter).  We  are  now 
left  with  choosing  p  in  (33)  to  minimize  the  variance  of  any  entry  1  <  b  j  <  cl- 

that  p  therefore  depends  on  ij.  This  parameter  p  is  in  practice  chosen  according  to  (15). 


3.4.2  A  second-order  model 

We  now  introduce  a  more  refined  approach,  that  not  only  takes  into  account  the  contributions 
from  single  defects  (through  see  (32))  but  also  contributions  from  pairs  of  defects.  To 

that  aim,  we  introduce 

=  (35) 
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_ I 

where  ,  defined  by  (30),  is  the  “marginal  eontribntion”  to  the  homogenized  matrix 

assoeiat<xl  to  the  eonfignration  with  tw’o  defects  located  in  A;  +  Q  and  I  +  Q.  In  view  of  (28), 
we  notice  that 


E 


l^Xi\i ,  l^k 


Hi 

2 


E  E  = 


,2-« 

2  ? 


kEX^  lEXp^ ,  l^k 


which  is  the  second  order  correction  in  the  expansion  (27).  When  77  is  small,  the  expectation 
ofX;„  +  2ir(  uj)  +  A2’^{uj)  is  a  good  approximation  of  the  expectation  of  A*  accurate 

up  to  an  error  of  the  order  of  rj^. 

In  a  way  similar  to  (33),  we  now  introduce  our  second-order  controlled  variable  as 

=  ^r,,v(^)  -  Pi  (^i'^^(^)  -  )  -  P2  (^2’^(w)  -  )  •  (36) 

We  have  introduced  two  deterministic  parameters  pi  and  p2,  which  need  not  be  equal.  For 
any  choice  of  these  parameters,  we  have  E  ~  [-^v  n]  • 

To  evaluate  (36),  we  first  have  to  precompute  the  deterministic  matrices 

-jk,N  _  -jO,N  -jk,l,N  _  -^,l-k,N 

^Idef  —  ^Idef  ^2def  —  ^2def  ' 


Computing  M  realizations  of  fherefore  amounts  to: 

•  offline  stage:  (i)  determine  by  solving  the  problem  (19)  (20)  on  Q  and  by  solving 

only  once  the  problem  (22)-(23)  on  Qv  (say  for  k  =  0);  (ii)  determine  by  solving 

|Ijv|  —  1  problems  (25)-(26)  on  (for  A;  =  0  and  /  G  J/v,  /  7^  0). 

•  online  stage:  solve  M  corrector  problems  (7)-(8)  on  (for  M  i.i.d.  realizations  of  A 

on  Qn)i  and  evaluate  according  to  (36). 


Questions  related  to  the  cost  of  evaluating  d.2 


are  discussed  below. 


Notice  that,  in  the  above  construction,  we  have  considered  as  reference  configuration  the 
defect-free  material,  i.e.  that  for  77  =  0.  Since,  in  the  regime  we  focus  on,  77  is  not  small, 
there  is  no  reason  to  favor  the  defect-free  configuration  (77  =  0)  rather  than  the  full  defect 
configuration  (77  =  1),  which  corresponds  to  the  periodic  matrix  Cper-  We  therefore  introduce 
(compare  with  (29)) 


— k  N 

^Idef  “  ~  ^per’ 


where  is  the  homogenized  matrix  corresponding  to  a  unique  defect  with  respect  to  the 

periodic  configuration  Cper  (compare  with  (21),  (22)  and  (23)): 


VpeR'',  CT,,,„p=W/' 

IQa'I  -Jqn 

where,  for  any  p,  the  corrector  is  a  solution  to 

-div  (Cf  (p-b  Vt;!  ^’^))  =0, 


Cl  (p  +  Vt;!-^’^)  , 

1,1, /c,v  jg  Q^-periodic, 


(37) 
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.  Likewise,  we  introduce  the  second  order  correction 

(38) 


where  Cf  =  Cper  -  Ifc+g^^C'per  - 
(compare  with  (30)): 

k,l,N  _ 

2def  —  '^2,k,l,N  "r  *-^per) 

where  Ch^  j^  is  defined  l)y  (37)  and  Clk,i,N  (conipare  with  (24),  (25)  and  (26)): 


Vp  £  R-*,  P  =  p^-r  /  e‘''  (p  + 

IVnI  ./q^ 


where,  for  any  p  €  K'^,  the  corrector  is  a  solution  to 

-div  {p  +  =  0,  is  Qyv-periodic, 

where  =  Cp^r  -  (u+g  +  L+q)  (c'per  -  >lper)-  As  in  (35),  we  introduce 

=  E  (l  -  B.»)  (l  -  B,») 

k^Xjsj  l^Xjsjyl^k 

where  C2dei  defined  by  (38).  Its  expectation  reads 


:=  E 


c:> 


.N 


=  5E  E  E  1(1  -  B")  (1  -  B,”)!  c‘ 

/cGXjv  l‘¥^k 


T^kylyN 

def 


=  5E  E 

fceXjv  ^GXtv,  l^k 


2  ■p5?Ai,/,A/' 

^  2  def  • 


We  eventually  introduce  the  controlled  variable  (compare  with  (36)) 


(39) 


(40) 


=  K.nM  -  PI 


VA^ 


-  PC  (vlj"(a,)  -  pMf )  -  P3  (c?’''(ai)  -  Vf)  .  (41) 


Consider  now  a  specific  entry  1  <  f,  J  <  d  of 
approach  consists  in  approximating  E  {A* j^) 


the  homogenized  matrix.  The  control  variate 
by  considering  a  Monte  Carlo  estimator  for 


3,r;  ') 

Pl,P2,P3/ 


The  deterministic  parameters  pi,  p2  and  are  ideally  chosen  to  minimize 


the  variance  of  P3(^))ij-  They  depend  on  ij  and  are  the  solution  of  the  following  3x3 

linear  system  (we  drop  the  subscript  i,j  for  conciseness): 


Var[A;'A]pi  +  CovIA;'’"",  A^’^]p2  +  Coy[A^{^,  C^'^]ps 
Cov[A''’"^,  Al^]px  +  Var[A^'^]p2  +  Cov[A^’^, 
CovlC^'"^,  +  Cov[a^'^,  +  Var[C2"A]p3 


Cov[yl*_^,Yl?’^] 

Cov[a*;v.A2^] 

^oy[A*N,C^’^] 


depending  on  the  covariances  between  the  entries  ij  of  A*  A^’^ ,  A^'^  and  C2'^ .  In  practice, 
these  covariances  are  approximated  by  empirical  estimators,  as  explained  in  Section  3.2. 
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In  practice,  computing  the  matrices  (and  likewise  is  rather  expensive  (because 

each  problem  is  set  on  the  large  domain  Qyv,  and  the  number  of  these  problems  increases  when 
N  increases).  It  is  therefore  useful  to  approximate  them  using  the  Reduced  Basis  strategy 
introduced  in  [16],  which  dramatically  decreases  the  computational  cost.  The  procedure  is 
essentially  as  follows.  We  first  solve  the  single  defect  problem  (22)  for  k  =  0,  and  solve  (25) 
for  a  limited  number  of  locations  of  the  defect  pairs,  say  k  —  0  and  I  close  to  k.  On  the  basis 
of  these  computations,  we  are  then  in  position  to  obtain  very  efficient  approximations  of  the 
matrices  ^2  def  ^  ^  I  7^  0.  Evaluating  (35)  is  thus  inexpensive.  Thus,  up  to  a 

limited  offline  cost  (i.e.  the  cost  for  solving  the  few  problems  (25)  that  we  have  to  consider), 
the  Monte  Carlo  empirical  estimator  and  the  Control  Variate  empirical  estiinal,or,  defined 
respectively  by 


rMC 


M 

—  V 


M 


{uj)  and 


rCV 

■'m 


M  ^  ‘ 


m=l 


/>1,P2,P3 


m=l 


share  the  same  cost. 


Remark  2  In  shanp  contrast  to  the.  first  order  control  variable,  the.  second  order  control  vari¬ 
able  not  only  depends  on  the  number  of  defects  in  the  materials,  i.e.  but  also  on 

keXN 

their  location.  The  specific  yeornetry  of  the  materials,  which  is  ignored  in  (34),  is  taken  into 
account  in  (41). 

Ill  the  secpiel,  we  di'iiioiistrate  iiuinerically  the  efficiency  of  the  approach  (see  e.g.  Figure  5 
below),  before  giving  some  theoretical  arguments  proving  that  variance  is  indeed  reduced. 


3.5  Numerical  experiments 

We  now  apply  the  methodology  described  above  to  some  two-dimensional  model  material  for 
which  the  field  A  is  of  the  form  (16)-(17)  (18)  (see  Fig.  3).  We  choose 

^per(2;)  =  «Id2  and  Cper(2;)  =  ,dld2, 

with  o;  =  3  and  /3  =  23  (similar  qualitative  conclusions  are  obtained  with  other  generic 
values).  All  variances  are  estimated  on  the  basis  of  M  =  100  independent  realizations.  All 
the  correctors  problems  have  been  solved  using  FreeFEM-f— 1-  [21],  on  a  mesh  of  size  h  =  0.2, 
using  PI  finite  elements.  We  prc'seiit  here  some  numerical  results,  and  refer  to  [6]  for  additional 
results. 

On  Fig.  4,  we  plot  as  a  function  of  G  (0, 1)  three  quantities: 

•  the  first  entry  of  the  matrix  E  [A*  (obtained  in  practice  by  an  expensive  Monte  Carlo 
estimation); 

•  the  weakly  stochastic  approximation  (27),  which  is  an  approximation  of  E  [A*  ^]  with 
an  error  of  the  order  of  0{rf)-, 

•  the  weakly  stochastic  approximation  obtained  in  the  regime  {1  —  p)  <C  1,  which  is  an 
approximation  of  E  [A*  with  an  error  of  the  order  of  0((1  —  p)^). 
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Figure  3:  A  typical  realization  of  the  checkerboard  test-case  with  r]  =  1/2. 


We  work  with  N  =  10,  and  the  following  observations  are  also  valid  for  larger  values  of  N. 
We  see  on  Fig.  4  that,  when  rj  <  0.4,  the  deterministic  expansion  (27)  is  a  very  accurate 

approximation  of  E  (A*  .  This  approximation  is  inexpensive  to  compute.  The  same 

observation  holds  in  the  regime  r]  >  0.7,  where  the  deterministic  expansion  around  rj  =  1 
provides  a  satisfying  approximation.  However,  we  note  that  none  of  the  two  weakly  stochastic 

expansions  are  accurate  when  0.4  <  77  <  0.7.  In  that  regime,  one  has  to  compute  E 
by  considering  several  realizations  of  (7)-(8).  In  that  regime,  considering  a  variance  reduction 
approach  is  useful. 


Supercell  size:  10x10 


Figure  4:  E  as  a  function  of  rj,  for  N  =  10.  Black  curves:  weakly  stochastic 

approximations.  Blue  curve:  Monte  Carlo  standard  estimator. 


In  the  regiiiK'  wc  have  identihed,  we  show  on  Fig.  5  the  ratios  of  variance 


7?r;,V  — 


Var(  [A;^] 
Var(Il) 


(42) 


where  D  is  <'ither  the  first-order  controlled  variable  Il^’''(a;)  defined  by  (33),  or  the  sc'eond- 
order  controlled  variable  D‘^p^p^{u})  dehned  by  (36),  or  the  controlled  variable 
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dofiliod  by  (41).  The  parameter  p  (resp.  (pi,P2)  and  (pi,p2ip3))  is  chosen  to  minimize  the 
empirical  variance  of  the  estimator  (see  Section  3.2). 

Remark  3  The  second-order  controlled  variable  defined  by  (36)  is  built  by  consid¬ 

ering  Aper  as  the  reference.  One  could  alternatively  build  a  second-order  controlled  variable 
considering  Cper  as  the  reference.  Numerical  results  obtained  with  such  a  controlled  variable 
are  similar  to  those  obtained  with  (results  not  shown). 


Supercell  size:  10x10 


a  I  ftn  QP  p/}  bd 


Figure  5:  Ratio  defiiK'd  by  (42)  as  a  function  of  g  {N  =  10).  Black  curve:  controlled 
variable  Dp'^{uj).  Red  curve:  controlled  variable  Blue  curve:  controlled  variable 

(uj). 

We  observe  on  Fig.  5  that,  for  g  =  the  approach  using  the  first-order  controlled 

variable  (33)  provides  a  variance  reduction  ratio  (42)  close  to  6.  This  gain  is  close  to  the  gain 
obtained  using  an  antithetic  variable  approach  (see  [13,  Table  2]).  In  contrast,  when  using 
the  controlled  variable  (41)  taking  into  account  first  order  and  second  order  corrections  with 
respect  to  both  the  cases  p  =  0  and  p  =  1,  we  obtain  a  gain  close  to  40. 

We  now  investigate  how  the  gain  depends  on  the  size  of  the  domain  Qj^.  To  that  aim,  we 
show  on  Table  1  the  ratio  (42)  as  a  function  of  N ,  for  p  =  1/2.  We  observe  that  the  gain  is 
essentially  independent  of  N . 


N  =  6 

N  =  10 

N  =  20 

II 

o 

N  =  50 

First  order 

7.57 

5.18 

6.55 

8.51 

7.34 

Second  order 

35.9 

41.8 

37.6 

35.6 

40.4 

Table  1:  Ratio  dehned  by  (42)  as  a  fuiK-tioii  of  N  (p  =  1/2).  First  order:  controlled 
variable  Second  order:  controlled  variable 

We  eventually  plot  on  Fig.  6  the  confideiice  intervals  obtained  for  the  Monte  Carlo  ap¬ 
proach  and  the  Control  Variate  approach  based  on  (41).  The  latter  confidence  interval  width 
is  dramatically  smaller  than  the  former. 
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Figure  6:  Estimation  of  as  a  function  of  N.  Blue:  standard  Monte-Carlo 

estimator.  Red:  Control  Variate  estimator  based  on  (41).  In  both  cases,  estimators  are  built 
using  M  =  100  i.i.d.  realizations. 

The  numerical  results  presented  above  have  been  exposed  in  details  in  [6].  They  clearly 
demonstrate,  from  a  numerical  perspective,  the  efficiency  of  the  approach. 


3.6  Theoretical  validation 

In  [6],  we  have  given  arguments  that  guarantee  the  efficiency  of  the  method  in  the  one  dimen¬ 
sional  case.  We  have  also  provided  there  elements  of  theoretical  analysis  in  higher  dimensional 
settings.  We  review  here  the  main  results  that  we  have  obtained. 


3.6.1  One-dimensional  case 

In  the  one-dimensional  case,  we  have  the  following  result: 

Proposition  1  Consider  the  model  (16)-(17)-(18)  in  the  one-dimensional  case.  Let  A* 
hr.  the  apparent  hornoqenized  matrix  defined  by  (7) -(8)  and  be  the  first-order-  controlled 
variable  defined  by  (33).  Then 


Var(A*  -  —  -\-0 


1 


and,  for  the  optimal  value  of  the  deterministic  parameter  p, 

min =  O  ■ 

Let  „  (^)  f>(-  the  seeond-order  controlled  'variable  defined  by  (41).  For  the  optimal  value 

of  the  deterministic  parameters  pi,  p2  and  pz,  we  have 

J £)3,r/  ^  -  t'  I 

\  PI, P2,P3  J 


min  Var( 

P1,P2,P3 


o 


iV3 
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Using  tlio  control  variate  approach  based  on  the  first-order  model  (resp.  second  order 
model),  the  variance  is  thus  improved  by  at  least  one  order  (resp.  two  orders)  in  terms  of  N. 
Note  in  particular  that,  in  the  above  results,  we  have  not  assumed  rj  to  be  small. 


3.6.2  Multi-dimensional  case 


In  the  multi-dimensional  case,  our  theoretical  results  cover  the  regime  1.  This  is  not  the 
regime  we  wish  to  address  with  our  approach.  Obtaining  theoretical  results  in  the  general 
case  is  still  an  open  question.  Note  that  the  numerical  experiments  of  Section  3.5  have 
been  perfored  in  the  regime  of  interest  when  weakly  stochastic  approximations  (based  on  the 
asymptotics  ??  -C  1)  do  not  yield  accurate  results. 


(■4U)„ 


can 


Consider  any  entry  ij  of  the  homogenized  matrix.  The  estimation  of  E 
be  done  by  a  Monte  Carlo  empirical  mean  on  {A*  i^  (see  Section  3.4.1)  or 

(see  Section  3.4.2). 


Proposition  2  For  any  entry  ij  of  the  homogenized  matrix,  we  have 


Yar 


=0(1,),  Var  (BjDy  =  Var  =  0{„’). 


We  thus  see  that,  by  using  the  first-order  (resp.  the  second  order)  surrogate  model,  the 
variance  is  improved  by  at  least  one  order  (resp.  two  orders)  in  terms  of  g. 


4  Homogenization  approach  for  the  numerical  simula¬ 
tion  of  periodic  microstructures  with  defects 

[Work  expanded  in  [12,  1,  2,  3[.[ 


The  homogenization  of  (deterministic)  non-periodic  systems  is  a  well-known  topic.  Al¬ 
though  well  explored  theoretically  by  many  authors,  it  has  been  less  investigated  from  the 
standpoint  of  numerical  approaches  (except  in  the  random  setting).  In  collaboration  with 
X.  Blanc  and  P.-L.  Lions,  C.  Le  Bris  has  introduced  a  possible  theory,  giving  rise  to  a  nu¬ 
merical  approach,  for  the  simulation  of  multiscale  non-periodic  systems.  The  theoretical 
considerations  are  based  on  earlier  works  by  the  same  authors  (derivation  of  an  algebra  of 
functions  appropriate  to  formalize  a  theory  of  homogenization).  The  numerical  endeavour 
that  is  on  the  horizon  is  completely  new. 

TIk'  general  approach  consists  in  approximating  at  the'  fine  scale  the  solution  to  an  elliptic 
ec[uatiou  with  <rseillatory  eoeffieieiit  wlic'ii  this  coefficient  c-oiisists  of  a  “nice”  function  which 
is,  in  some  sense  to  be  made  precise,  perturbed.  A  typical  example  is  that  of  a  periodic 
material  (where  the  period  is  small  with  respect  to  the  size  of  the  sample,  so  that  there 
exists  a  separation  of  scales  between  the  microstructure  and  the  macrostructure),  with  a 
superimposed  defect  in  one  (or  ”a  few”)  periodic  cells.  The  equation  of  interest  is 


— div 


+  Bdef 


=  /  ill 


(43) 
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where  the  perturbation  is  localized  (in  the  sense  that  B^ei  vanishes  at  infinity),  while 
Aper  encodes  the  perfect,  periodic  medium.  From  the  macroscopic  standpoint,  the  overall 
behavior  might,  or  not,  be  the  same  as  that  of  a  perfect  periodic  material.  Anyhow,  close  to 
the  defect  (close  to  the  ’’support”  of  B^ef),  thc'  response  is  of  course  veny  diffenait  from  that 
of  the  perfect  material. 

We  aim  at  d('V('loping  numerical  approaches  abb'  to  efficicaitly  and  accurately  capture  this 
local  difference.  In  terms  of  material  science',  this  is  of  course'  e)f  paramount  importance: 
material  failure  indeed  often  occurs  close  to  defects. 

The  approach  is  based  on  the  determination  of  a  local  profile,  solution  to  an  equation 
similar  to  the  corrector  equation  (2)  in  classical  homogenization.  In  the  case  (43)  with  Bjef  G 
L^(IR^)  (an  assumption  which  somehow  formalizes  the  fact  that  Sjef  vanishes  at  infinity),  this 
equation  is 

-div  [(Aper  +  Bdef)  [P  +  VtCp)]  =  0  ill  , 

where  p  is  any  fixt'd  vector  in  The  well-posedness  of  that  equation,  in  various  functional 
settings  depending  upon  the  nature  of  the  perturbation,  is  a  theoretical  issue  that  has  now 
been  investigated  and  solved,  for  several  prototypical  situations;  local  perturbation,  two 
diffen'iit  pi'riodic  stnicturi's  separati'd  by  a  common  int  erface,  . . . 

The  theoretical  results  obtained  to  date  are  being  collected  in  the  series  of  publications  [12, 
1,2].  The  work  [12]  also  contains  preliminary  numerical  simulations,  in  the  case  5def  G  //^(M'^). 
The  review  article  [3]  presents  the  various  approaches  within  a  more  general  perspective. 

The  purpose  is  now  to  consider  cases  closer  to  actual  local  defects  in  materials  (such  as 
dislocations)  and  also  to  put  all  this  in  action  as  a  general  numerical  approach  on  cases  of 
actual  practical  relevance. 


5  A  parameter  identification  problem  in  stochastic  ho¬ 
mogenization 


[Work  expanded  in  [7].] 

Consider  the  highly  oscillatory  equation  (4),  where  the  matrix  A  is  assumed  to  be  ran¬ 
dom  and  stationary.  We  have  seen  in  Section  2.2  that  this  problem  admits  the  homogenized 
limit  (3).  Random  homogenization  theory  actually  provides  formulas  to  compute  the  ho¬ 
mogenized  matrix  A*,  see  (5)  (6).  We  thus  have  at  our  disposal  a  procedure  to  compute 
macroscopic  quantities  if  we  know  the  microscopic  quantities,  and  to  solve  the  so-called  for¬ 
ward  problem.  However,  in  practice,  given  a  heterogeneous  matc'rials,  it  is  a  diffifuilt  question 
to  decide  on  the  law  of  the  microscopic  physical  properties,  i.e.  on  the  probability  law  of 
A{x,u)).  On  the  other  hand,  macroscopic  quantities  are  more  easily  accessible.  It  is  thus  of 
interest  to  consider  the  inverse  problem,  and  try  to  extract  some  information  on  the  properties 
of  the  materials  at  the  microscopic  scale  on  the  basis  of  macroscopic  quantities. 

Of  course,  homogenization  is  an  avc'raging  process,  which  filters  out  many  features  of 
t  he  microscopic  coefficic'uts.  There'  is  thus  no  hope'  to  recover  a  full  information  about  the 
microstructure  (in  our  case,  the  probability  distribution  of  A{x,uj))  from  the  only  knowledge 
of  macroscopic  quantities  such  as  A*.  We  adopt  here  a  more  restricted  objective.  We  assume 
a  functional  form  for  the  distribution  of  tlu'  microscopic  field  (uauu'ly,  a  Weibull  distribution. 
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see  Section  5.1  below)  and  aim  at  recovering  the  parameters  of  that  microscopic  law  of  the 
basis  of  macroscopic  quantities.  This  question  belongs  to  the  wide  family  of  inverse  problems. 
However,  the  niutiscale  context  brings  in  soiiu'  specificities.  We  refeu  to  [24]  for  a  review 
article  on  inverse  problems  in  such  a  multiscale  context. 

5.1  The  microscopic  model 

We  conskh'r  here  a  problem  written  at  the  fine  scale  as  a  finite  difference  problem  (rather 
than  as  the  partial  differential  etiuation  (4)).  The  motivation  for  this  choice  is  related  to 
the  particular  application  (transport  phenomena  in  porous  media,  modelled  as  a  lattice  of 
channels)  we  have  in  mind.  Howc'ver,  our  approach  is  not  specific  to  discrete  elliptic  equa¬ 
tions.  It  could  also  be  applied  for  problems  modelled  by  continuous  elliptic  partial  differential 
eciuations  (PDEs)  with  random,  highly  oscillatory  coefficients,  such  as  (4). 

For  the  sake  of  completeness,  we  detail  here  the  homogenization  result  in  the  setting  of 
discrete  elliptic  equations,  which  is  in  essence  the  same  as  that  recalled  in  Section  2.2. 

Let  P  be  a  bounded  domain  of  and  /  €  C^{T>).  Let  A  be  the  random  stationary  matrix 
field  given  by 

Vx  G  A{x,  uj)  =  diag(^ai{x,  uj), ...  ,ad{x,uj)j ,  (44) 

where  the  conductances  {aj(x,  •)}i<i<(i  form  an  i.i.d.  sequence  of  random  variables.  Let 
«£  be  the  unique  solution  to 

V*[v4(x/£,a;)VeMe(x,  cn)]  — /(x)  inPrieZ*^,  Ue(x,  cn)  =  0  in  (R*^  \  P)  fl  eZ*^,  (45) 

where  the  discrete  gradient  (resp.  the  discrete  divergence  V*)  is  defined,  for  any  function 
V  defiiu'd  on  tin*  lattice'  eZ'^,  by 


(Vet;)(x) 


v{x  -f  eei)  —  u(x)^ 

yu(x  -f  ecd)  -  v{x)  J 


and 


-(V»(x) 


Vi{x  -  eej) 
e 


The  equation  (45)  plays  here  the  role,  in  the  current  setting,  of  (4). 

When  e  goes  to  0,  converges  to  some  homogenized  function  u*,  solution  to  the  (contin¬ 
uous)  partial  differential  ecpiation  (3),  where  the  homogenized  matrix  A*  is  given  as  follows: 


Vp  G  A* p  =  E[A{x,-){p +  'Vwp{x,-))] 


(46) 


where,  for  any  p  G  R*^,  the  corrector  Wp  is  the  unique  (up  to  the  addition  of  a  constant) 
solution  to 


/ 


V* 


yl(-,a;)(p  +  V'U)p(-,a;)) 


=  0  in  Z'^, 


I'Vwp  is  stationary, 

Vx  G  Z'^,  E[VrCp(x,  •)]  =  0. 


(47) 
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The  expressions  (46)-(47)  play  here  the  role,  in  the  current  setting,  of  the  expressions  (5)-(6). 

As  in  the  continuous  case  described  in  Section  2.2,  the  corrector  problem  (47)  is  untractable 
in  practice,  since  it  is  posed  in  the  entire  lattice  Approximations  are  therefore  in  order. 
The  standard  procedure  amounts  to  considering  finite  boxes  as  in  Section  3.1.  For  any  N,  we 
denote  by  Qat  the  finite'  box  {0, . . . ,  N}^.  The  truncated  corrector  {-,00)  is  defined  on  Qat 
as  the  solution  to 

(  V7*  A/.  ,  I  _  n  /O.. 

(48) 


(49) 


f-v* 

■.4(-,w)(p  +  V<(-,u.))' 

l<(-. 

a;)  is  Qvperiodic. 

In  turn,  the  homogenized  matrix  A*  is  approximated  by  (hTned  by 

1 


Vp  G 


Wp  =  IQ 


N\ 


A{x,u){p  +  VWp{x,uj)). 

xEQj^ 


In  the  sequel,  we  assume  that  the  conductances  {o.i{x,uj)}^^2d  i<i<d  (entering  the  micro¬ 
scopic  fh'ld  A,  see  (44))  form  an  i.i.d.  sequence  of  random  variables  that  are  distributed 
according  to  the  Weibull  law  of  parameter  We  recall  that  such  random  variables  are 

positive,  with  a  probability  density  that  reads  (see  Figure  7) 


Vr  >  0,  =  X  (D  (“(^/^)^)  • 

In  practice,  a  \Wibull  distribution  is  generated  as  follows.  Let  u{cu)  be  a  random  variable 
uniformly  distributed  in  [0, 1].  Then 


fc-i 


i/k 


a(uj)  —  A  —  ln(l  —  u(uj)) 
is  distributed  according  to  the  Weibull  law  of  oarameter  (A.  k). 


(50) 


Weibull  (fistributions 


Figure  7;  Examples  of  \Wibull  distributions. 


5.2  Forward  and  inverse  problems 

The  forward  (direct)  problem  can  be  phrased  as  follows:  given  the  parameters  (A,  k)  of  the 
Weibull  law  of  the  microscopic  field  A(-,w),  compute 
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•  the  macroscopic  (homogenized)  permeability  E  , 


•  its  relative  variance 


VarR[A;v]  = 


Var[A^] 


The  inverse  problem  we  consider  here  is  the  following.  Given  the  quantities  E  and 
VarR[74^],  we  wish  to  hnd  parameters  (A,  A")  consistent  with  these  observations.  The  reason 
why  we  choose  E  [v4^]  and  VarR[j4^]  as  input  parameters  is  discussed  below. 

Ideally,  the  inverse  problem  can  be  formulated  as  follows.  Consider  the  functional 

=  (EbM  _  AC  (y-mw)] 

V  ^obs  /  V  ^bs 

where 

^obs  ^  [-^7v(^obs7  ^obs)] 

is  the  observed  macroscopic  permeability,  while 


^^bs  (^obs  ?  ^obs)] 

is  the  observed  relative  variance.  We  consider  the  problem 


inf  Fat  (A,  k) 

X,k 

and  we  obviously  see  that  (Aobs,  A;obs)  is  a  minirnizer  of  the  above  problem.  Furthermore,  we 
have  the  following  result  (see  [7]): 

Lemma  1  Consider  the  one- dimensional  setting  d  =  1,  and  set 

Fooi^,k)=  lim  FN{\,k). 

N—^OO 

Then  Foo(A,A;)  has  a  unique  minirnizer,  which  is  (Aobs,  A^obs)- 

The  above  result  motivates  our  choice  of  the  expectation  and  the  relative  variance  as  input 
observations.  It  shows  that,  although  homogenization  is  an  averaging  process  which  filters 
out  many  features  of  the  microscopic  coefficients,  our  two  macroscopic  input  observations  are 
enough  to  characterize  the  two  parameters  of  the  microscopic  probability  distribution. 

In  practice,  the  expectations  in  are  approximated  by  empirical  estimators: 

1 

E  [/l;v(A,  A;)]  k',  oj)  ;=  —  A,  k) 

m=l 

and  likewise  for  VarR[/l)^(A,  A-)],  which  is  approximated  by  IAv,m(A,  k\(jj).  A  practical  formu¬ 
lation  is  hence  to  consider  the  optimization  problem 


inf  Fv,a/(A,  A;w), 

A,fc 


(51) 


where 


F,  u) 


/4* 

^obs 


+ 


f  Fiv,M(A,  A;  u) 

V  Kbs 


where  the  observed  values  are  given  by  =  A^jv,^(Aobs,  Aobsi^*^)  and  likewise  for  Kbs- 
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5.3  Algorithm  and  numerical  results 

To  address  the  optimization  problem  (51),  recall  that,  in  view  of  (44)  and  (50),  the  microscopic 
hold  can  bo  writt('ii  as 

4l(’)  ^)  =  fc,  {ufc(a;)}^g2‘i  j ) 

where  Uk{Lo)  are  i.i.d.  random  variables  uniformly  distributed  in  (0, 1),  and  where  ^  is  a 
deterministic  function,  the  derivatives  of  which  (with  respect  to  A  and  k)  are  easy  to  compute. 

By  computing  the  derivatives  of  (48)  and  (49)  with  respect  to  A  and  k,  one  can  easily 
oouiputo  the  hrst  and  second  derivatives  of  k;Lo).  It  is  thus  possible  to  use  a  Newton 

algorithm  to  solve  (51). 

On  Figure  8,  we  show  the  numerical  results  that  we  have  obtained  on  a  two-dimensional 
problem,  with  the  parameters  A  =  10  and  M  =  30.  We  proceed  as  follows.  We  pick  some 
target  values  Aobs  and  fcobs,  pick  some  hxed  cJ  G  fl,  and  solve  the  forward  problem  for  these 
values,  thereby  generating  synthetic  observed  values  jy^(Aobs)  ^’obs)^)  and  Kbs.  We 

next  repeat  the  following  procedure; 

1.  we  generate  microstructures  A{x,u)  at  the  hiu'  scale  whidi  are  indepeiidc'iit  of  those 
used  to  compute  the  observed  values  and  Vobs- 

2.  we  solve  the  optimization  problem  (51)  using  a  Newton  algorithm,  starting  from  an 
initial  guess  10%  off  the  reference  value  (Aobsi^obs)- 

3.  we  thereby  identify  some  optimal  parameters 

We  repeat  the  above  procedure  many  times  in  order  to  obtain  several  i.i.d.  realizations  of  the 
optimal  parameters  (Aopt(n;),  k^pt{<^)),  from  which  we  build  the  histograms  shown  on  Figure  8. 
We  see  that,  despite  the  limited  values  of  N  and  M,  we  obtain  a  meaningful  estimation  of 
the  exact  parameters  (Aobs  )  ^obs  )  • 


,  .  r  kovti'^)  ■”  ^obs  T->-  1  .  1  •  J.  C  '^opt(^)  Aobs 

Figure  8:  Left;  histogram  of  ^  ^  - .  Right;  histogram  of - — - . 

Remark  4  Our  <i'i)p‘i'0(ich  is  not  sj)ccific  to  Weihull  Uiws.  It  can  be  uscdfoT'  other  distribution 
laws  with  parnmeters  0.  Whut  we  need  is  that  the  ‘mndoni,  fi(dd  A(^x,uj)  used  at  the  microscopic 

scale  can  be  written  as  .  . 

A{x,u})  =  IF(u{x,ijo),0] 
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where  u(x,  u)  is  a  field  of  random  variables  that  arc  uniformly  distributed  and  J-  smoothly 
depends  on  the  parameters  6  (see  (52)  in  our  particular  case).  Computing  the  derivatives  of 
the.  microscopic  random,  field  A{x,ui)  (and  next  of  the  macroscopic,  homogenized  quantities) 
with  respect  to  6  is  then  easy.  Our  motivation  for  choosing  Weibull  laws  comes  from  physical 
reasons:  based  on  experimental  results,  it  appears  to  be  a  reasonable  choice. 

6  Conclusions 

The  series  of  works  presented  in  this  report  has  contributed  to  provide  more  efficient  nnmerical 
methods  for  the  simulation  of  random  heterogeneous  materials. 

First,  in  some  previous  works  (see  [13]  and  subsequent  works),  funded  by  an  earlier  ONR 
grant,  we  have  demonstrated  the  feasibility  of  using  variance  reduction  techniques  in  the  con¬ 
text  of  stochastic  homogenization.  We  used  there  the  antithetic  variable  approach.  We  have 
shown  in  Section  3  that  it  is  as  well  possible  to  use  other  variance  reduction  approaches,  such 
as  control  variate  approaches,  that  involve  the  construction  of  reduced  models  in  this  particu¬ 
lar  context.  This  technique  is  less  generic  than  the  antithetic  variable  approach,  but  provides 
better  results.  Ongoing  efforts  are  focused  on  the  d(^v(dopmeut  of  yet  another  technique, 
based  on  the  a  priori  selection  of  “representative”  microstructures  [4].  This  selection  step  is 
inexpensive  in  comparison  to  the  overall  computational  cost.  The  corrector  problem  (8)  is 
next  solved  only  for  these  “better”  microstructures. 

Of  course,  considering  other  problems  than  the  linear  equation  (4)  (such  as  nonlinear  prob¬ 
lems,  . . . )  is  also  of  interest.  In  [5],  we  have  considered  a  nonlinear,  convex  homogenization 
problem,  and  we  have  shown  that  the  technique  of  antithetic  variables,  considered  in  [13]  for 
linear  problems,  can  also  be  used  in  that  nonlinear  context,  with  similar  results. 

Second,  we  have  started  to  develop  homogenization  approaches  for  the  numerical  sim¬ 
ulation  of  periodic  microstructures  with  defects  (see  Section  4).  This  theory  gives  rise  to 
completely  new  numerical  approaches.  Numerical  results  have  already  been  obtained  for 
some  classes  of  defects.  These  questions  have  not  been  explored  in  the  previous  ONR  grant. 
Our  purpose  now  is  to  consider  cases  closer  to  actual  local  defects  in  materials. 
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